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0. TielemanQ A. Lazarides, and C. Morais Smith 
Institute for Theoretical Physics, Utrecht University, 
Leuvenlaan 4, 3584 CE Utrecht, The Netherlands 
(Dated: September 16, 2010) 

We present the theoretical mean-field zero-temperature phase diagram of a Bose-Einstein con- 
densate (BEC) with dipolar interactions loaded into an optical lattice with a staggered flux. Apart 
from uniform superfluid, checkerboard supersolid and striped supersolid phases, we identify several 
supersolid phases with staggered vortices, which can be seen as combinations of supersolid phases 
found in earlier work on dipolar BECs and a staggered-vortex phase found for bosons in optical 
lattices with staggered flux. By allowing for different phases and densities on each of the four sites 
of the elementary plaquette, more complex phase patterns are found. 
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I. INTRODUCTION 

During the last decade, the field of ultracold atomic 
gases in optical lattices has undergone tremendous 
growth PQ, having witnessed amongst others the reali- 
sation of the supernuid-Mott insulator (SF-MI) transi- 
tion [2j [3], the construction of spin-dependent lattices 
[4], and the generation of synthetic gauge fields [5]. In 
the earlier years, substantial effort was aimed at real- 
ising and understanding the Bose- and Fermi-Hubbard 
models, which feature only on-site interactions. Recently, 
there has been much progress in bringing longer-ranged 
interactions into the system by using dipolar atoms or 
molecules [6 . An example of experimental progress 
in Bose-Einstein condensation of dipolar atoms can be 
found in Ref. [7 . Another approach is the creation of 
dipolar molecules; see e.g. Ref. [8]. Theoretical studies 
have indicated numerous interesting possibilities, includ- 
ing, but not limited to, supersolidity [6, 9-11 . 

Supersolidity is commonly defined as the simultaneous 
presence of diagonal and off-diagonal long-range order in 
the system [12], a prominent candidate for experimen- 
tal realisation being solid 4 He p~3j[14]. Other candidates 
have been suggested in the domain of ultracold atomic 
gases, such as rapidly rotating Fermi-Fermi mixtures [T5] 
and dipolar bosonic gases in optical lattices [10]. In an 
optical lattice, one clearly has diagonal long-range or- 
der, since the density at the minima of the lattice po- 
tential is higher than at the maxima; however, this type 
of long-range order is imposed externally. To preserve 
the analogy to bulk supersolids, where the diagonal or- 
der is spontaneously present in the system, we define a 
supersolid in an optical lattice as a phase with both long- 
range off-diagonal and diagonal order, where the diago- 
nal order breaks the translational symmetry of the lat- 
tice [16] . Recently, dipolar atoms or molecules in optical 
lattices have been predicted to feature such supersolid- 
ity. Refs. [TTJ [17] present analytical and numerical analy- 
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ses of dipolar atoms in square lattices with only nearest- 
neighbour (NN) hopping. Other examples include square 
lattices with NN and next-nearest-neighbour hopping [18] 
and triangular lattices [T6] , 

Uniform magnetic fields for ultracold atomic gases have 
been mimicked by applying rotation [I5j [19] and by gen- 
erating gauge fields using position-dependent optical cou- 
pling between internal states of the atoms [20 , 21 . Anal- 
ogous to superconductors in magnetic fields, these sys- 
tems exhibit vortices. A staggered gauge field for neutral 
atoms has also been proposed [22] , leading to a staggered- 
vortex superfluid phase [23]. In this paper, we find that a 
dipolar bosonic gas subjected to a staggered gauge field 
exhibits a supersolid phase which features vortices. In 
contrast to Ref. [15] , we study the gas in a lattice and do 
not have a rotating trap. 

We analyse the interplay between NN interactions 
and an artificial staggered magnetic field in a system of 
bosons in a two-dimensional square optical lattice. In 
order to perform this analysis, we generalise and com- 
bine the methods used in Refs. [10], [23], and [24], to 
allow for the description of phases with a higher degree 
of broken symmetry than discussed in those three refer- 
ences. We present a phase diagram containing combina- 
tions of the uniform and staggered-vortex phases found 
in Ref. [23] and the supersolid phases found in Ref. [10] , 
as well as a region where two phases coexist and the sys- 
tem will phase separate. We find that several continuous 
and discontinuous phase transitions between different su- 
perfluid and supersolid phases can be driven in two ways: 
by changing the NN interaction strength or by changing 
the applied flux. Apart from the presence or absence 
of density modulations, we discuss the existence of an- 
other type of structure in the system, which arises when 
the many-body wavefunction exhibits phase differences 
between neighbouring lattice sites. The vortices studied 
in Refs. 23 and [24 are a realisation of a non-trivial, 
although relatively simple phase structure. 

This paper is structured as follows: In section [TTJ we 
introduce the system and briefly discuss its constituent 
components. In section [TTT| we present the methods used 
to determine the phase diagram, which is displayed and 
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discussed in sections [IV] and El In section ED we show 



experimental signatures of the phases found. Section [VII 
concludes the paper by summarizing and discussing the 
results. 



II. THE SYSTEM 

We consider a system of dipolar bosons in a two- 
dimensional square optical lattice [10] with staggered flux 
[24] . Below, we briefly explain the consequences of a stag- 
gered flux (s ectio n II A ) and dipolar interactions in a lat- 
tice (section |lIB ) for the phases found in the system. We 
work at T = and only consider Bose-condensed phases, 
since we are interested in the combined effects of a stag- 
gered flux and anisotropic NN interaction in a superfluid. 
The Hamiltonian we investigate is 



H — influx + i^on-site + ^dip 

= " J E (e i( - 1)V 4& r+e! +H.c.) 



U 



^ n r (n r - 1) 



(1) 



+ ^2^2 {y x n r n r+aGl +V y n r n r+(je ^j. 

reAcr=±l 

Here, J represents the hopping strength between neigh- 
boring sites; is the phase picked up at each jump, which 
is related to the magnitude of the flux through a pla- 
quette; U the on-site interaction strength; and V x and 
V y the anisotropic NN interaction strengths. The lattice 
is represented as two interspersed square sublattices, A 
and B. The operators a Y (a\) and b r / are destruc- 
tion (creation) operators acting on sites in the sublat- 
tices A and B, respectively; note that there is only one 
type of particle being created and destroyed. The oper- 
ator n r is the number operator for site r, irrespective of 
the sublattice in which it is located. The lattice vectors 
ej, I G {1,2,3,4} are defined by ei = — e3 = e x and 
e 2 = -e 4 = e y . 



A. Staggered flux 

The term Hr ux breaks the symmetry between the sub- 
lattices A and B, as can be seen in Fig. 1(a). It can 
be represented as a synthetic magnetic field which alter- 
nates in sign between neighboring plaquettes. For the 
details of the derivation of this hopping term, we refer 
the reader to Refs. [23] and [22] . The phase diagram of 
bosons with on-site interactions in such a lattice is pre- 
sented in Ref. [23]; we reproduce it in Fig. 1(b). The 
main conclusions from Ref. [23] which will be important 
for this paper are the following: For strong on-site inter- 
actions (large U), a Mot t- insulating phase is found. By 
reducing J7, a second-order phase transition into a su- 
perfluid phase occurs at some critical value of U/J; the 
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FIG. 1: (color online), (a) The division into sublattices A 
and B caused by the staggered flux, (b) The phase diagram 
for bosons with on-site interactions in an optical lattice with 
staggered flux as found by the authors of Ref. [23]. Legend: 
SF = conventional superfluid; SVSF = staggered- vortex su- 
perfluid; SSSF = staggered-sign superfluid. 



value of (U/J) c depends on as shown in Fig. 1(b). For 
|0| < 7r/4, i.e. small flux, the zero- momentum superfluid 
phase is unaltered; for 7r/4 < |0| < 37r/4, the system 
features a staggered- vortex superfluid phase, where the 
vortex cores are located at the centers of the plaquettes 
and the sign of the vorticity alternates between plaque- 
ttes. This phase comes about due to the development 
of a second minimum in the single-particle spectrum, at 
a finite momentum, which becomes the global minimum 
if 7r/4 < |0| < 37r/4. Condensation in this minimum 
introduces phase differences of tt/2 between the lattice 
points, in such a pattern that a particle tunnelling around 
a plaquette picks up a phase of ±27r, depending on the 
direction of tunnelling and the plaquette. Note the peri- 
odicity in the 0-dependence of (U/J) c . We find that the 
same periodic pattern emerges in the V x -(j)- diagrams that 
we present in sections [IV| and \V\ respectively, in spite of 
the fact that we are studying different phase transitions. 
We will confine ourselves to the weakly interacting limit, 
since our aim here is to study the interplay between the 
supersolidity found in Ref. [10] and the staggered- vortex 
patterns found in Ref. [23] . 



B. Dipolar interaction 

The interaction energy between two polarized dipoles 
is given by 



Vdd(r) = d 2 g dd 



1 -3cos 2 C 



where £ is the angle between the polarisation axis and the 
displacement vector r. Loading the dipoles into a deep 
lattice and approximating the dipolar interaction by cut- 
ting it off at NN distance, the only displacement vectors 
that we have to consider are e x and e2, the lattice vectors 
in the x- and ^-directions. The NN dipolar interaction 
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strengths in the two relevant directions are given by 

_2 1 — 3 sin 2 d cos 2 (p 
V x =d g dd = , 

T 3 



a° 

1 — 3 sin 2 $ sin 2 ip 



where a is the lattice spacing, $ is the inclination, and 
(p is the azimuthal angle. At cp = 7r/4, the interaction 
strength is isotropic and can be varied continuously from 
repulsive (# = 0) to attractive (# = 7r/2), being zero at 
d = sin -1 y/2/3. By varying the azimuthal angle, we can 
tune the ratio between V x and V y to any desired value. 
We note that tuning the NN interactions to zero will 
make the next-nearest-neighbor interactions more rele- 
vant; however, in this paper, we focus on the strong NN 
interaction regime. Since the relevant quantities for the 
purposes of our analysis are V x /U and V y /U : we can 
cover the complete V x /U-V y /U -plane by tuning J7, and 
cp. Thus, we simply have two independent dimensionless 
parameters, V x /U and V y /U. The Hamiltonian for the 
dipolar interaction thus takes the form 



(2) 



r,a=±l 



In Ref. [10], the phase diagram of dipolar bosons in a 
square optical lattice is presented, which we reproduce in 
Fig. 2 and briefly discuss here. At mean-field-level, the 
authors identify three types of superfluid phases: a con- 
ventional superfluid one (SF), for weak NN interactions; 
one with a density modulation in a checkerboard pattern 
(checkerboard supersolid - CSS), for strong enough re- 
pulsive V x ~ V y ; and one with a striped pattern (striped 
supersolid - SSS), where only one of the NN interac- 
tion parameters dominates. Both the checkerboard and 
striped superfluids have long-range diagonal as well as 
off-diagonal order, i.e. they break both gauge and trans- 
lational symmetry. These properties justify the epithet 
supersolid, as discussed in the introduction. 
The checkerboard and striped phases are intuitively easy 
to understand. Consider the case where V x = V y : 
here, the NN interaction energy is reduced by arranging 
the atoms in a checkerboard-modulated density pattern, 
since there are fewer pairs of nearest neighbors in such 
a configuration. Similarly, if V y is negative (or positive 
but small) and V x is positive and of the order of Z7, the 
NN interaction energy is reduced by a striped configura- 
tion. As a final remark, taking longer-range interactions 
into account does affect the phase diagram to some ex- 
tent, but the basic structure of checkerboard and striped 
phases remains intact [TO] . 
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FIG. 2: (color online) The phase diagram for dipolar bosons 
in a square optical lattice. The black and blue lines were also 
found by the authors of Ref. 10 . Legend: SF = homogeneous 
superfluid, CSS = checkerboard supersolid, SSS = striped 
supersolid. SSSl (SSS2) has the stripes in the y-(x-) direction. 
The solid black lines are second-order phase boundaries, the 
dashed black lines first-order phase transitions, and the dotted 
blue (grey) lines represent the existence of metastable states, 
which are labelled in blue (grey) and between brackets. The 
red (grey) dashed lines are our finding, and represent the SSS- 
CSS phase boundary for the case of strong flux, cj) = 7r/2; see 
the discussion in section [V] 



Bogolyubov approximation to describe the condensate, 
which is justified if J ^> U,V X: V y [10 ; we calculate the 
ground state energies and excitation spectra of all the 
phases we identify. The excitation spectrum is required 
to check the dynamical stability of the phases we find: 
since we are dealing with bosons, a negative or complex 
excitation spectrum implies that the ground state above 
which the spectrum is calculated is not the real ground 
state of the system. 

We will be working in the limit of high filling factor v 
and strong hopping (or weak interaction), such that we 
can consider every lattice point to contain a condensate, 
and still have a small J / (yU), where v is the filling factor. 
The ratio J/ (vU) is required to be small for the density 
modulations to appear; if the energy gain from wavefunc- 
tion overlaps between neighboring sites is too large, the 
system will ignore the NN interactions and simply remain 
in a superfluid state with homogeneous density. 



Two-sublattice formalism 



III. THE METHOD: MEAN-FIELD 

Below, we explain how we obtain the phase diagram 

and [V] We use the 



that we present in sections IV 



As described in section II A[ the staggered flux divides 
the square lattice into two interpenetrating square sub- 
lattices, necessitating a two-sublattice description of the 
system, which is developed below. In the next subsection, 
we introduce another subdivision, into four sublattices, 
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in order to allow for more complex density and phase 
distributions around the elementary plaquette. 

We will perform our calculations in momentum space, 
where the grand-canonical version of the Hamiltonian 
given in Eq. ([!]) reads 

H — fiN = ^2 ( e k«k^k + h.c. - /iaj^k - nb\.b\^j 



kGBZi 



E 



+k 2 -k 3 -k 4 ,0 



ki,...,k 4 
€BZi 



u 

2N« 



(4x4 



k 2 a k 3 «k4 



2F(k 4 



K^k 4 ) 



a kl &k 2 a k 3 &k 4 



(3) 



where 



a: a_ k 



hi b. 



and cjk, Ak, 7k, and (k are functions to be calculated for 
each specific phase we describe with this formalism. We 
diagonalise this quadratic Hamiltonian by solving 



(M k -fi k [A k ,4]) =0, 



for ^k, which then is the excitation spectrum. 



with 

e k = -4 J 



/ \ / H~ fey 
cos((/>)cos' 



k^jQ ky 

a I cos I — -a 



■ . / i . i k x H~ ky \ . / k x ky 
- 2 sin (0 J sin I — -a | sm( — -a 



and 



V(k) = V x cos(e x • k) + V y cos(e y • k). 



N is the particle number operator for the entire system, 
H is the chemical potential, N s is the number of sites per 
sublattice, i.e. half the total number of sites, and BZi is 
the first Brillouin zone. Now, we apply the Bogolyubov 
approximation: we replace — )► (5k, c (a) + ^k and —> 
^k,c (b) + ^k, where c is the condensation momentum, 
and we treat and 6k as small fluctuations relative to 
the average occupations (a) and (b). We then require 
that the terms linear in the fluctuations vanish. This 
requirement yields the values for the chemical potential 
and, if present, the density modulation; for the details, 
see sections IIVI and [Vj The terms of order zero in the 
fluctuations represent the ground state energy, 



£ =2Re(e c (a) (b)*) + U(\ (a) | 4 - 
+ 2n0)|(a)| 2 |(6)| 2 /iV s . 



(b) | 4 )/2iV s 



(4) 



The second-order terms can be diagonalised to give the 
excitation spectrum. They can be represented in matrix- 
form as 



keBZi 



Ak 7k Ck " 
K ^k c k 7 k 
7k Ck w k A k 
- Ck 7k A£ w k 



AlM k A k , 

keBZi 



(5) 



B. Four-sublattice description 



The method presented above is in fact a simple com- 
bination of the approaches used in Refs. [10] and [23] . 
It works as long as either the phase distribution around 
the elementary plaquette is trivial (i.e. all sites have the 
same phase) or the density modulation is absent. If we 
allow V x 7^ V y , striped or otherwise asymmetric phases 
may occur, in which case the four sites of the elementary 
plaquette could all have different densities, or the phase 
drops could be distributed unevenly along the plaquette. 
In such phases, the two-sublattice formalism does not 
hold, since the density modulation and the phase pat- 
tern will influence each other, as we show below. One 
interesting phenomenon to be investigated here is the 
competition between NN interactions favouring stripes 
on the one hand, and a staggered flux on the other hand, 
since the staggered flux is associated with c heckerboard- 
subdivision of the lattice (see section II A Fig. 1). To 



investigate such phenomena, we need a description of 
the system which allows the condensate wavefunction to 
be different at all four sites of the elementary plaquette. 
Such a description involves four different sublattices: we 
have to split each of the sublattices A and B into two 
new ones, such that A = SL10SL3 and B =SL2®SL4: 



a b a b 1212 
baba^434 3' 



The resulting equations for the chemical potential and 
condensate wavefunction cannot be solved analytically, 
and have to be solved numerically instead. In this repre- 
sentation, the momentum-space Hamiltonian becomes 
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H = ^ ( e k a i,k a 2,k + ( 6 k)* a i,k a 4,k + £k4,k a 4,k + (4)* a 3,k a 2,k) + h.c. 



kGBZi 



" V ^2 5Z a lk a *,k+ ^ ^k 1 +k 2 -k 3 -k 4 ,0 5Z a lk 1 a i,k 2 a *,k3 a i,k 4 
kGBZi i=l ki,...,k 4 ^ 2=1 

A/" 

- cos[(k 4 - k 2 ) • e y ] (4, kl 4,k 2 a i,k 3 a 4,k 4 + 4,^4,^ a 2,k 3 a3,k 4 ) 



■ COs[(k 4 - k 2 ) • e x ] (4,k 1 4,k 2 a l,k 3 «2,k 4 + 4,k 1 4,k 2 a 3,k 3 «4,k 4 ) 



(6) 



r 



where 



conditions 



u/y 



2 J cos(k x / y a)e 



Note that this four-sublattice formalism can in princi- 
ple be used to describe all the phases that we will en- 
counter in the system, including the ones that can be 
analysed within the two-sublattice formalism; for exam- 
ple, the checkerboard- modulated phase without vortices 
has (ai) — (as) = /3 (a 2 ) = /3 (a 4 ), where j3 is a measure 
for the strenght of the modulation. 



1. Mean field 

In the four-sublattice description, the real space unit 
cell is twice as large as in the two-sublattice description. 
As a consequence, the first Brillouin zone is only half 
the size, and the number of bands in the excitation spec- 
trum is doubled. The minimum in the corner of the first 
Brillouin zone found in Ref. [23] is mapped to the cen- 
ter of the new first Brillouin zone. Hence, the minimum 
of the single-particle spectrum in the four-sublattice de- 
scription is always at k = 0, and we use the mean- field 
ansatz a^k — » #k,o (aj) + a J? k- We obtain four equations 
for the chemical potential, of the form 

^ = ~ 2J ( fll ) - + ^'(ai)| 2 (7) 

+ 2V x is\(a 2 )\ 2 + 2V y is\(( H }\ 2 

(we omit the other three for brevity; they can easily be 
deduced from the Hamiltonian). These four jij can be 
interpreted as the chemical potentials of the four sublat- 
tices, which are then thought of as macroscopic systems 
with an exchange mechanism (the hopping terms). The 
condition that all four chemical potentials are equal rep- 
resents the equilibrium condition in this picture. Since 
the expression in Eq. ^\ is in principle complex, we have 
to allow for complex {aj ) in order to be able to make 
the imaginary parts of the jij vanish. Representing (aj) 
as r^e ld \ the requirement that all \ij are real yields the 



r 2 sin a\ = r 4 sin a 4 
rs sin a 2 = n sin a\ 
r 4 sin as — r 2 sin a 2 
r\ sin a 4 = r% sin a^, 



(8) 



where aj = Oj+i — 0j + (j). The four equations (|8| can 
be reduced to three without loss of generality. A fourth 
equation comes from the requirement that the mean-field 
wavefunction is always single- valued. To satisfy this re- 
quirement, the phase picked up when hopping around a 
plaquette has to be an integer multiple of 2tt; thus, 



ai — A<j) = 27m, 



(9) 



where n determines the vorticity pattern of the system. 
Apart from Eqs. (J8|, we also have the real parts of the 
chemical potentials, which have to be equal to each other. 
They take the form 

2 J 

/ii = (r 2 cos a\ + r 4 cos a 4) + U vr\ 



ri 
2V x vr. 



(10) 



2V y vri 



and similar expressions for /i 2 , ^3 and /i 4 . Finally, there 
is the normalisation condition, 



ri+ri + ri + ri = 4, (11) 
to a total of eight equations in eight un- 



bringing us 

knowns. Of course, more than one unique solution may 
still exist. We have to find all unique solutions, since each 
represents a different phase of the system, and we want to 
compare all phases for stability and ground state energy. 
Hence, we solve the eight equations Eqs. (|8|-(11) numer- 
ically, in such a way that we find all solutions. Having 
found the (%•), the ground state energy is given by 

4JRe J2 ^ M* <%+i> + ^ £ I M I" 



£n = 



2V, 



f (I (ox) | 



(« 2 ) 
(a 4 ) 



2N 

j 

<a 3 > | 2 | <a 4 ) | 2 ) 
(«3>| 2 |(a 2 )| 2 ) 



(12) 
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where j is to be taken modulo 4. 



2. Excitation spectrum 

As before, finding all solutions that represent equilib- 
rium situations is not enough: we also have to investigate 
their excitation spectra to assess their respective dynam- 
ical stabilities. In order to derive the excitation spectra, 
we collect all terms of second order in the fluctuations, 
and obtain, quite generally, 



ff « =2 S[ w k4k^k + ^*4k a ^ 



kGBZi i,j 



(13) 



«i,k%,-k 



where 



^k 



a l,k a l,-k <4,k a 2,-k <4 )k a 3,-k «4 5 k «4,-k 



In our specific case, cu u and X u do not depend on k: they 
are given by 



22 

^k 

,33 



/i 
/i 



^ = 217| (ax) | 2 + 2^| (a 2 ) | 2 + 2V y \ <a 4 > | 
:2/7|(a 2 )| 2 + 2^|(a 1 )| 2 + 2^|(a 3 )| 
a£ 3 = 2tT| (a 3 ) | 2 + 2V X \ (a 4 ) | 2 + 2^| (a 2 ) | 2 - /i (14a) 
^=2U\(a 4 ) 

xi 



■2V x \(a 3 )\ 2 + 2V y \(a 1 )\ 2 -ii 



:/7((a,)*) 2 . 

The hopping terms have the following coefficients: 



In most regions of parameter space, there is only one 
dynamically stable phase. However, there are some re- 
gions where two phases are dynamically stable; here, we 
compare the ground state energies (which are equal to 
the free energies, since we are working at T = 0) to see 
which phase is favoured. 

The phase diagram we find is, in principle, three- 
dimensional, since we have the three parameters (/>, V x 
and V y . In addition to the V^/ £7- "^/[/-diagram at (j) = 
from Ref. [10], we present three cross sections which to- 
gether form a representative sample of the results: the 
14-0-diagram at V y = V x , the V^-^-diagram at V y = 0, 
and the V x /U-V y /U -diagram at <p = tt/2. 



IV. QUANTUM PHASES: SYMMETRIC CASE 

In this section, we discuss the symmetric case, where 
the NN interaction is equally strong in both directions. 
This can be achieved by polarising the dipoles perpen- 
dicular to the plane, or at an angle of 7r/4 relative to the 
in-plane lattice vectors. By tuning the inclination, the ra- 
tio V x / y /U can be tuned without changing U. However, 
if this technique is employed, the next-nearest-neighbor 
interactions will not be isotropic; this is only the case if 
the polarisation axis is perpendicular to the plane. As a 
consequence, the description presented here will be more 
accurate if the dipoles are polarised perpendicularly to 
the plane. 



, ,23 



^ = u 3 r = 



.41 



,13 



^k 



: K )* 



e£ + 2V X cos(k x a) (ai) (a 2 )* 
e k + 2V y cos(k y a) (a 2 ) (a 3 )* 
e£ + 2V X cos(k x a) (a 3 ) (a 4 )* 
e k + 2V v cos(k y a) (a 4 ) (ai)* 

0. 



(14b) 



A. Weak flux: no vortices 

For small values of the flux, < \(j)\ < 7r/4, the system 
does not feature any vortices in the ground state. It may 
still exhibit density modulations, since these are caused 
by the NN interactions. 



Lastly, the off-diagonal mixing terms read 



A k 


- A k 


= 2V X 


cos(k x a) (ai)* (« 2 )* 


\23 
A k 


a 32 
- A k 


= 2Vy 


cos(k y a) (a 2 y (a 3 )* 


A k 


- A 43 

- A k 


= 2V X 


cos(k x a) (as)* (a 4 )* 


A k 


- A 14 

— A k 


= 2V y 


cos(k y a) (a±)* (ai)* 


A k 


- A k 


= 




A k 


- A k 


= 0. 





1. Superfluid: homogeneous density 

For homogeneous or checkerboard-modulated density 
distributions, we can use the two-sublattice formalism, 
and replace ak ^ 4 o (a) + &k and — » 5^ o (b) + b^. 

(14c) 

If the density is homogeneous, we have a conventional 
superfluid (SF) and | (a) \ = | (b) \ = y/N p /2, with N p 
being the total number of particles in the system. Setting 
the term linear in the fluctuations zero yields 



As before, we diagonalise this Hamiltonian by numeri- 
cally solving 



for ^£ sl , which then gives the excitation spectrum. 



fJ.SF = eo + v(U + 2V x +2V y ). 

Note that eo = — 4 J cos 0, and hence as <p — > 0, we recover 
the result found in Ref. [10]: /i S F = -47 + v(U + 2V X + 
2V y ). In order to determine where this superfluid phase is 
dynamically stable, we consider the excitation spectrum 
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given by Eq. (|5| . The matrix elements of M^ F are 

w k = 2v(U + V X + V y ) - fisF 

\ k =vU 

7k = e k + Ck 

C k = 2z/y(k). 

As (j) — >• 0, the excitation spect rum should reduce 
to nl F = V / e k [e k + 2^(/7 + 2V(k))], with e k = 2J[2 + 
cos(k x a) + cos(fc 2/ a)], as found in Ref. [10]. However, 
in order to make this comparison properly, we have 
to map the two-band excitation spectrum in the Bril- 
louin zone defined by k x ± k y G [— n/a, ir/a] derived 
above, to the single-band spectrum in the Brillouin zone 
k x / y G [— 7r/a,7r/a] given in Ref. [ID]. This mapping is 
described in e. g. Ref. [25 : it is the mapping from the ex- 
tended to the reduced zone scheme. After applying this 
mapping, the two spectra are seen to be identical. In 
addition, as V x — >• and V y —> 0, the spectrum reduces 
to the one found in Ref. [251. 



Now, we can now write a and (3 in terms of the parame- 
ters J, (/>, U, V x and V y : 

1 2Jcosd> 



a = 



2 1/(2^ + 2^-17)' 
1 2Jcosd> 



2 1/(2^ + 2^ - f7) 

With this result, we can calculate the chemical potential, 
and find ficss = 2vU . For the checkerboard-modulated 
case, the matrix elements of M^ ss in Eq. ([5| are given 

by 

a£ /6 = 2Uv{l ± 2^3) + 2V(0)i/(l T 2^) - Mess, 
A^ /6 =C/z/(l + 2v^), 
7k =£k + Ck, 

Ck=2z/y(k). 

As in the homogeneous case, if we take the limit 0—^0, 
we recover the results from Ref. [10] : density modulation, 
chemical potential, and excitation spectrum. Note that 
for < \(j)\ < 7r/4, the density modulation is affected by 
the flux, even though no vortices appear (see Eq. (15)). 



2. Checkerboard supersolid 



B. Strong flux: staggered- vortex phase 



To describe a phase with a checkerboard-modulated 
density (checkerboard supersolid, CSS), we follow the 
scheme used in Ref. [10] . We assume that the population 
density of sublattice A is different from that of sublat- 
tice S; hence, (a) ^ (b). Instead of working with (a) 
and (b) directly, we go to a center-of-mass and relative 
representation, and write 



(a)=JN p /2(V^ + 



(b}=^jN p /2(V^-VP). 

Clearly, if we send j3 to zero, the density modulation van- 
ishes. Replacing a k —> £ k ,o (o) + and 6 k — > 5k, o (b) + 6 k 
in Eq. Q and requiring the terms of first order in the 
fluctuations to vanish yields 

AU _^o V / Q-V? 
^ v ^/a - 

where v = N p /N s is the average number of atoms per 
site (filling factor) and we have used the fact that that 
| (a) | 2 + | (b) | 2 = N p and hence a + /3 = 1. Setting 
= I^B yields a condition on a and /3, which can be 
solved for the difference a — /3; the result is 



17 + 2V(0) + 2y/a/3[U - 2V(0)], 
U + 2V(0) - 2Vc^[/7 - 2V(0)], 



a-/3 = 



4 J cos ( 



1/(2^ + 2^-17)" 



(15) 



As was found in Ref. [23 , under the influence of a 
strong flux, 7r/4 < \<p\ < 37r/4, the system goes to a 
staggered-vortex superfluid (SVSF) phase. This corre- 
sponds to condensation in the single-particle state with 
momentum (±7r/a, 0) or (0,±7r/a), i.e. in the corners 
of the first Brillouin zone [26 . To describe this region 
of parameter space, we need the four-sublattice formal- 
ism presented in section [Til B[ Although the description 
is different, the ansatz is still informed by earlier find- 
ings: we expect a combination of the staggered-vortex 
pattern from Ref. [23 and, for appropriately strong NN 
interactions, the checkerboard- modulated density from 
Ref. [10 . For a homogeneous density distribution, the 
ansatz is quite simple: 

where we have defined the phase of the mean-field wave- 
function at sublattice SLi to be zero. For the dis- 
cussion of the staggered-vortex checkerboard supersolid 
(SVCSS), there is a general point that is worth noting: 
in cases where the density modulation is invariant under 
exchange of the two lattice vectors, the vortices do not 
interfere with the density modulation. This can be seen 
from Eq. (J8|: as long as T\ = and V2 = ^4, the condi- 
tions on the phase drops do not depend on the wavefunc- 
tion amplitudes. Since the checkerboard pattern has this 
symmetry, we can try to guess the modulation strength 
from the two-sublattice formalism. The ansatz would be: 



«k -> yj JVp/2(5 k ,7r {y/a + y^) + a k 
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Performing the same analysis as for the weak-flux case, 
we find that the density modulation strength is given by 



a — j3 



4 J sin ( 



v{2V x + 2V y - U) 



(16) 



Hitherto, the two formalisms work equally well. However, 
the excitation spectrum can only be calculated in the 
four-sublattice formalism; the mean-field values for the 
four sublattices are 



<ai> = sjN v /2{yfa 
<a 2 ) = - i v / iVp/2(VS 
<«3> = - 



N p /2(^ + 
<a 4 > =i v / iV2(^-^). 



By inserting these values for (aj) into Eqs. (12) and (14), 
we find the corresponding ground state energy and exci- 
tation spectrum. 



C. Phase diagram 

Now we have all the information required to deter- 
mine the cross section of the phase diagram along the 
line V x = V y . In order to obtain the phase boundaries, we 
calculate the parameter values where the density modula- 
tion vanishes, where the excitation spectra become unsta- 
ble, and where the ground state energies of two phases are 
equal. These three phenomena happen simultaneously 
at the SF-CSS and SVSF-SVCSS phase boundaries (see 
Fig. 3). Since the relevant order parameter for this phase 
transition is the density modulation, which vanishes con- 
tinously at the border, this is a second-order phase tran- 
sition. At the SF-SVSF and CSS-SVCSS boundaries, 
the ground state energies of the uniform and staggered- 
vortex phases become equal. The relevant order param- 
eter for these transitions is the (staggered) number of 
vortices per plaquette, which jumps from zero to unity; 
hence, these are first-order phase transitions. 

Fig. 3 shows the V x = V y cross section of the phase 
diagram. The combination of flux and NN interactions 
leads to a phase boundary with the same shape as ob- 
served in Ref. [23] in the MI-SF phase diagram. The 
phase transition from the homogeneous-density phase to 
the checkerboard phase is continuous, independently of 
the flux; the phase transition from the uniform to the 
staggered- vortex phase is discontinuous, independently 
of the density modulation. 

In conclusion, this cross-section of the phase diagram 
shows a relatively straightforward superposition of the 
phases found in Refs. [10 and [23]. Beyond the results 
found in those references, we note a few aspects: The 
symmetry in the line (j) = 7r/4 goes further than the shape 
of the phase boundary, since the strength of the density 



modulation is also symmetric. The origin of this symme- 
try can be found in Eqs. (15) and (16): the cosine from 
the uniform SF and CSS phases is replaced by a sine in 
the staggered-vortex phases. The SF-CSS transition can 
also be induced by changing 0, provided V x and V y are 
in the right range. The critical value ofV x /U shifts with 
J '/vU ", but the phase boundaries retain their periodicity 
in 0, within the approximation used here. 

?r/2 r 



(f>n/4 



SVSF 


svcss 


SF 


1 CSS 


1/4 


1/2 



Vy and 



FIG. 3: Cross section of the phase diagram for V a 
J/(yU) — 0.1. Dashed (solid) lines are first (second) order 
phase transitions. Legend: SVSF = homogeneous staggered- 
vortex superfluid; SVCSS = staggered-vortex checkerboard 
supersolid. 



V. QUANTUM PHASES: ASYMMETRIC CASE 



In Eqs. 



and (10), we see the interplay between 
As noted in see- 



the phase and density distributions, 
tion |IVB[ both a homogeneous density and a checker- 
board pattern are possible within the ansatz employed 
in Ref. [23], in which all phase drops are assumed to be 
equal. If, however, r\ ^ r% and/or r2 ^ r±, as in the 
case of stripes or four different densities on the four sites 
of one plaquette, the phase distribution is influenced by 
the density modulation, and we have to allow for phase 
drops taking other values than integer multiples of tt/2. 
We indeed find such solutions, in the striped phase for 
< <t> < 7r/2. The mean- field wavefunction values (aj) 



are determined by numerically solving Eqs. (pj)-(ll). 



A. Phase diagram cross section I: V y = 

We obtain the phase diagram in the same manner as in 
section [IV| by calculating where the density modulations 
vanish, where the excitation spectra develop instabilities, 
and by comparing ground state energies. Along the line 
V y = 0, we find a variety of phases (see Fig. 4). The 
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FIG. 4: Cross section of the phase diagram for V y = and 
J/OC/) = 0.1. Legend: SSS = striped supersolid; SVSSS = 
staggered- vortex striped supersolid; FR = forbidden region, 
phase separation. Dotted blue (grey) lines represent the ex- 
istence of metastable states, which are labelled in blue (grey) 
and between brackets. 



superposition of the CSS phase from Ref. [10 with the 
staggered- vortex phase from Ref. [23 is there, as well 
as the uniform and staggered- vortex superfluids without 
density modulation. The striped supersolid (SSS) and 
staggered- vortex striped supersolid (SVSSS) phases are 
not simple combinations of earlier found phases, however: 
they break the symmetry within the sublattices A and S, 
i.e. between sublattices SLi and SL3, and SL2 and SL4, 
respectively. The density modulates in a striped pattern, 
but the phase drops are different on all four edges of 
the elementary plaquette. Note again the high degree of 
symmetry in the line <p = tt/4: the dynamical stability 
diagram is completely symmetric, as well as the density 
modulation strength, also for the striped phase. Apart 
from the phase distribution, which cannot be symmet- 
ric, it is only the ground state energy difference between 
the checkerboard and striped phases that is different be- 
tween the two regions \(j)\ < 7r/4 and \(j)\ > 7r/4. In the 
staggered- vortex region, the checkerboard phase has a 
lower ground state energy than the striped phase, which 
can be understood as a consequence of the matching be- 
tween the sublattice divisions associated with the flux 
and the density modulation. The flux breaks the symme- 
try between sublattices A and B (for details see Ref. [23 J, 
thus introducing a checkerboard pattern, which competes 
with the striped pattern introduced by the NN interac- 
tions. 

Also note the shape of the second-order phase bound- 
ary, in this case between the homogeneous and striped 
phases. It shows the same pattern as the boundary be- 
tween the homogeneous and checkerboard phases (see 
Fig. 3), and the SF-MI boundary in the absence of NN 
interactions (see Fig. 1). This shape can be understood 



from the effect of the flux on the hopping energy in the 
ground state: it is modified by cose/) in the uniform su- 
perfluid, and sin0 in the staggered- vortex superfluid, re- 
sulting in a minimum at <ft = 7r/4. Since it is the hopping 
term in the Hamiltonian that favors the superfluidity in 
the SF-MI phase transition, and the homogeneity in the 
SF-CSS and SF-SSS transitions, the reduction in hopping 
energy makes phases which break the phase coherence or 
homogeneous density distribution more favorable. 

There is one region, close to (j) = 7r/4, where we do 
not find any dynamically stable phases (FR - forbidden 
region). This result is a consequence of our approach, 
which assumes the existence of a well-defined chemical 
potential for the whole system, and hence a uniform 
macroscopic density distribution, since we are not work- 
ing in a trap. If we drop the assumption of a well-defined 
chemical potential, we effectively allow the system to sep- 
arate into parts with different densities. Shifting the den- 
sity changes the parameter J ' jvJJ which changes the crit- 
ical values of V x /U in such a way that any point in the 
forbidden region can be made to lie within the two clos- 
est neighboring stable regions. Hence, if the parameters 
are tuned to lie in the forbidden region, our calculations 
predict that the system will phase separate. Note that 
there are actually two different forbidden regions, since 
the phases the system will separate into are different for 
|0| < 7r/4 and |0| > 7r/4. 



B. Phase diagram cross section II: (j) — tt/2 

Finally, a comment on the two cross sections of con- 
stant (j) (see Fig. 2). We show the cases where = 
and (j) = 7r/2, i.e. the centers of the uniform and the 
staggered-vortex phases. As mentioned above, the phase 
diagram has a high degree of symmetry in the line 
(j) = 7r/4, the only differences being the phase distribu- 
tion and the ground state energies of the striped and 
checkerboard phases. In Fig. 2, we see that the CSS-SSS 
phase boundary shifts, but nothing else. Intermediate 
cross sections would also reveal the forbidden region dis- 
cussed above, and the disappearance of the striped phase 
near <\> = 7r/4. 



VI. EXPERIMENTAL SIGNATURES 

A good starting point for experimental detection of 
the various phases discussed in this paper is the momen- 
tum distribution n(k), since most of the phases have a 
unique momentum distribution, as will be discussed be- 
low. Experimentally, the momentum distribution is ac- 
cessible through the technique of time-of-flight measure- 
ment [1 , which converts the momentum distribution into 
a spatial one by suddenly turning off the lattice and al- 
lowing the cloud to expand ballistically. n(k) is given by 
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(a) 



1 







FIG. 5: The momentum distributions for the (a) SF, (b) CSS, 
and (c) SSS phases. 



FIG. 6: The momentum distributions for the (a) SVSF, (b) 
SVSSS at cj) = 7tt/20, and (c) SVSSS at = tt/2 phases. The 
SVCSS phase has exactly the same signature as the SVSF 
phase. 



n(k) = \w(k) 



E^ kP £ 



ik-(r 



«v0 



Here, P runs over all plaquettes, i.e. P = 2a(n x ,n y ), 
with integer n x / y , and r M and r v run over the four sites 
of the plaquette. The last factor is the one that distin- 
guishes between the phases, since all configurations are 
invariant under translation by P. To calculate it, we fol- 
low Ref. [23] and approximate the TV-particle state on the 
lattice by a coherent state with on average N particles. 
Now the calculation becomes very simple, since in this 
approximation, 



for which we have explicit results (see sections TTTJV). 
Using the same parameters as in Ref. [23] , we obtain the 
signatures shown in Figs. 5 and 6. In Fig. 5, we see that 
supersolidity manifests itself in the momentum distribu- 
tion by replacing the peaks from the homogeneous SF by 
a smaller peak and 'satellite peaks' displaced by the char- 
acteristic vectors of the density modulation. Fig. 6 shows 
that the same replacement takes place in the staggered- 
vortex phases. Unfortunately, this implies that the SVSF 
and SVCSS phases look exactly the same, since the peaks 
of the SVSF momentum distribution are displaced from 
each other by exactly the characteristic vectors of the 
CSS density modulation. In Fig. 6, we see another cu- 
riosity: the two striped phases, with stripes in the x- and 
^/-directions, are indistinguishable for <j> = tt/2, where the 
phase drops along the plaquette are equal. For different 
values of </>, the two striped phases are distinguishable. 
Also, as discussed in sections |III| and [V] the phase dis- 
tribution is asymmetric for <\> £ {0,7r/2} in the striped 
phase. This asymmetry reflects the interplay between 
NN interactions and applied flux, and is a continuous 
function of 6. 



VII. DISCUSSION & CONCLUSIONS 

In this paper, we have analysed the interplay between 
NN interactions and a synthetic staggered magnetic field 
in a system of bosons in a two-dimensional square op- 
tical lattice. We have used the Bogolyubov approxima- 
tion to obtain the theoretical mean-field phase diagram of 
the system. The equilibrium condition that traditionally 
gives the value of the chemical potential was replaced by 
a set of conditions that give the density and phase mod- 
ulations between the lattice sites, as well as the chemical 
potential. The excitation spectrum allowed us to deter- 
mine the dynamical stability of the phases encountered 
in the system. 

Our analysis resulted in a rich phase diagram featur- 
ing various superfluid and super solid phases. Apart from 
the conventional and staggered-vortex superfluids and 
the checkerboard and striped supersolids found before, 
the system turns out to feature phases which combine 
a staggered-vortex phase configuration and supersolid- 
ity. Where the density modulation is invariant under 
exchange of the lattice vectors, it does not influence the 
phase configuration; in case of an asymmetric density 
modulation, the phase drops around the plaquette are 
also distributed asymmetrically. Lastly, we have identi- 
fied a forbidden region, where the system cannot form a 
stable state with the same average density in all areas. In 
this region of parameter space, our calculations predict 
that the system will phase separate. 

We observe that even a rather crude approximation 
to the dipolar interaction, where the long tail is cut 
off beyond the NN range, combined with the staggered 
flux, leads to a very rich phase diagram. Apart from 
the expected 'superposition' of density modulations and 
a staggered-vortex pattern, another layer of structure 
emerges: the phase drops do not have to be distributed 
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homogeneously along the elementary plaquette. 

We also see that many of the phase transitions can 
be driven by tuning either the NN interaction strengths 
or the flux, which is a consequence of the fact that the 
flux modifies the hopping energy. Thus, it affects both 
the vorticity, which is a discrete variable, and the density 
modulation, which is a continuous variable. This is yet 
another example of the interesting physics that comes 
with the possibility of generating an artificial staggered 
magnetic field in an optical lattice. 

Potentially interesting questions that were beyond the 
scope of this paper include taking into account the effects 
of the long tail of the dipolar interaction and the effects 
of finite temperatures on the supersolid phases. Since 
our work was exploratory in nature, we have not been 
able to address these problems, but they would certainly 
be relevant for experimental tests of the predicted phases. 
Another point left unaddressed here is the full periodicity 



in <p of the phase diagram. Fig. 1 shows the 27r-periodic 
nature of the system, while we have only considered val- 
ues of cj) between and tt/2 in this paper. Note that 
since the two staggered-vortex phases are identical, the 
system is symmetric under <p — >> —cj). Thus, we have stud- 
ied half of the phase diagram's entire period. The other 
half is expected to feature phases and phenomena similar 
to those already discussed here, as can be deduced from 
the nature of the density and phase distributions found. 
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